ICASE 


DISCRETE-TIME  MARKOVIAN  STOCHASTIC 
PETRI  NETS 


Gianfranco  Ciardo 


19950329  038 


Contract  No.  NAS  1-19480 
February  1995 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


Operated  by  Universities  Space  Research  Association 


Discrete-time  Markovian  stochastic  Petri  nets 

Gianfranco  Ciardo 

Department  of  Computer  Science 
College  of  William  and  Mary 
Williamsburg,  VA  23187-8795 
ciardo@cs.wm.edu 

Abstract 

We  revisit  and  extend  the  original  definition  of  discrete-time  stochastic  Petri  nets, 
by  allowing  the  firing  times  to  have  a  “defective  discrete  phase  distribution”.  We  show 
that  this  formalism  still  corresponds  to  an  underlying  discrete-time  Markov  chain.  The 
structure  of  the  state  for  this  process  describes  both  the  marking  of  the  Petri  net  and 
the  phase  of  the  firing  time  for  of  each  transition,  resulting  in  a  large  state  space.  We 
then  modify  the  well-known  power  method  to  perform  a  transient  analysis  even  when 
the  state  space  is  infinite,  subject  to  the  condition  that  only  a  finite  number  of  states 
can  be  reached  in  a  finite  amount  of  time.  Since  the  memory  requirements  might  stiU 
be  excessive,  we  suggest  a  bounding  technique  based  on  truncation. 
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NASA  Contract  No.  NAST19480  while  in  residence  at  the  Institute  for  Computer  Applications  in  Science 
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1  Introduction 

In  the -past  decade,  stochastic  Petri  nets  (SPNs)  have  received  much  attention  from  re¬ 
searchers  in  the  performance  and  reliability  arena  and  have  been  extensively  applied  to 
the  performance  and  reliability  modeling  of  computer,  communication,  manufacturing,  and 
aerospace  systems  [4,  5,  7,  10,  23].  While  there  is  agreement  on  the  appropriateness  of  SPNs 
as  a  description  formalism  for  a  large  class  of  systems,  two  radicalN  different  solution  ap¬ 
proaches  are  commonly  employed:  simulation  and  state-space-based  analysis.  Simulation 
allows  to  associate  general  distributions  to  the  duration  of  activities  (SPN  transitions),  but 
it  requires  multiple  runs  to  obtain  meaningful  statistics.  This  problem  is  particularly  acute 
in  reliability  studies,  where  many  runs  might  be  required  to  obtain  tight  confidence  intervals. 
With  simulation,  the  state  of  the  SPN  is  composed  of  the  marking,  describing  the  structural 
state  of  the  SPN,  and  the  remaining  firing  times,  describing  how  long  each  transition  in  the 
SPN  must  still  remain  enabled  before  it  can  fire.  The  simulated  time  9  is  advanced  by  firing 
the  transition  with  the  smallest  remaining  firing  time. 

State-space-based  analysis  has  been  mostly  applied  to  SPNs  whose  underlying  process  is 
a  continuous-time  Markov  chain  (CTMC),  that  is,  to  SPNs  with  exponentially  distributed 
firing  times  [3,  12,  25,  26].  Except  for  numerical  truncation  and  roundoff,  exact  results  are 
obtained,  but  the  approach  has  two  limitations:  the  number  of  SPN  markings  increases 
combinatorially,  rendering  unfeasible  the  solution  of  large  models,  and  generally-distributed 
activities  must  be  modeled  using  “phase-type  (PH)  expansion”  [15].  PH  distributions  can 
approximate  any  distribution  arbitrarily  well,  but  it  is  difficult  to  exploit  this  fact  in  practice 
because  the  expansion  exacerbates  the  state-space  size  problem. 

Discrete  distributions  for  the  timing  of  SPNs  have  received  less  attention.  This  is  un¬ 
fortunate,  since  deterministic  distributions  (constants)  are  often  needed  to  model  low-level 
phenomena  in  both  hardware  and  software,  and  the  geometric  distribution  is  the  discrete 
equivalent  of  the  exponential  distribution  and  can  approximate  it  arbitrarily  well  as  the  size 
of  the  step  decreases.  Furthermore,  there  is  evidence  supporting  the  use  of  deterministic 
instead  of  exponential  distributions  when  modeling  parallel  programs  [1]. 

If  all  the  firing  distributions  are  geometric  with  the  same  step,  the  underlying  process  is 
a  discrete-time  Markov  chain  (DTMC)  [2-5].  Such  SPNs  can  model  synchronous  behavior, 
as  well  as  the  main  aspect  of  asynchronous  systems:  the  uncertainty  about  the  ordering 
of  quasi-simultaneous  events.  A  DTMC  is  described  by  a  square  one-step  state  transition 
probability  matrix  H  and  an  initial  state  probability  vector  7rf°k  The  state  probability  vector 
at  step  n  can  be  obtained  with  the  iteration  [power  method):  =  tt  [n  This  result 

was  extended  in  [11]  to  include  immediate  transitions,  which  fire  in  zero  time,  and  geometric 
firing  distributions  with  steps  multiple  of  a  basic  unit  step,  possibly  with  parameter  equal 
one.  that  is,  constants.  [29]  restates  these  results  in  more  detail,  and  uses  the  concept  of 
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weight  to  break  the  ties,  following  [3]  and,  more  closely,  [13].  Generalized  Timed  Petri  Nets 
(GTPN)  have  also  been  proposed  [19],  where  the  steps  of  the  geometric  firing  times  for 
each  transition  can  be  arbitrary,  unrelated,  real  numbers.  A  DTMC  can  be  obtained  by 
embedding,  but  the  analysis  is  restricted  to  steady-state  behavior  and  the  state  space  of  the 
DTMC  can  be  infinite  even  when  the  underlying  untimed  PN  has  a  finite  reachability  set. 
Analogous  considerations  hold  for  D-timed  PNs  [30]. 

We  generalize  and  formalize  the  results  in  [13]  and  show  how,  using  phase-expansion,  a 
DTMC  can  be  obtained  even  if  the  firing  time  distributions  are  not  geometric,  as  long  as 
firings  can  occur  only  at  some  multiple  of  a  unit  step.  The  state  can  then  be  described  by 
the  marking  plus  the  phase  of  each  transition.  This  extends  the  class  of  SPNs  that  can  be 
solved  analytically,  but  two  limitations  still  exist:  the  existence  of  a  basic  step  and  the  size  of 
the  state  space.  By  using  a  fine  step,  arbitrary  steps  can  be  approximated,  but  this  increases 
the  state  space. 

Approaches  to  solve  models  with  a  large  state  space  have  been  proposed  for  both  steady- 
state  and  transient  analysis.  [6]  considers  the  reliability  study  of  a  SPN  with  exponentially 
distributed  firing  times,  under  the  condition  that  the  reachability  graph  is  acyclic.  The 
underlying  CTMC  is  then  acyclic  as  well,  and  a  state  can  be  discareded  as  soon  as  the 
transitions  out  of  it  have  been  explored,  resulting  in  an  algorithm  offering  large  savings 
in  memory  and  computations  with  respect  to  traditional  numerical  approaches.  However, 
acyclic  state  spaces  arise  only  in  special  cases,  such  as  reliability  models  of  non-repairable 
systems. 

For  transient  analysis  of  a  general  CTMC,  Jensen’s  method  [21],  also  called  uniformization 
[17,  27],  is  widely  adopted.  [18]  outlines  a  dynamic  implementation  of  the  algorithm,  where 
the  state  space  is  explored  as  the  computation  of  the  transient  probability  vector  proceeds, 
not  in  advance,  as  normally  done.  This  allows  to  obtain  a  transient  solution  even  if  the  state 
space  is  infinite,  provided  that  the  transition  rates  have  an  upper  bound. 

If  the  CTMC  contains  widely  different  rates,  the  number  of  matrix- vector  multiplications 
required  by  uniformization  can  be  excessive.  Proposals  to  alleviate  this  problem  are  selective 
randomization  [24]  and  adaptive  uniformization  [28],  both  based  on  the  idea  of  allowing 
different  uniformization  rates,  according  to  the  set  of  states  that  can  be  reached  at  each 
step.  The  latter,  in  addition,  can  be  used  with  infinite  state  spaces  even  if  the  rates  have  no 
upper  bound.  However,  the  method  can  incur  a  substantial  overhead,  and  it  appears  that 
an  adaptive  step  is  advantageous  only  in  special  cases  or  for  short  time  horizons. 
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In  Sections  2,  3,  and  4  we  define  the  underlying  untimed  PN  model,  the  class  of  DDP 
distributions  used  for  the  temporization  of  a  PN,  and  the  resulting  DDP-SPN  formalism, 
respectively.  Section  5  discusses  the  numerical  solution  of  a  DDP-SPN,  by  building  and 
solving  its  underlying  stochastic  process,  a  DTMC.  Section  6,  examines  approaches  to  cope 
with  large  state  spaces. 


2  The  PN  formalism 

We  recall  the  (extended)  PN  formalism  used  in  [12, 14].  A  PN  is  a  tuple  (^P,  T,  D~  y,g, 

where: 

•  P  is  a  finite  set  of  places,  which  can  contain  tokens.  A  marking  g  £  defines 

the  number  of  tokens  in  each  place  p  ^  P,  indicated  by  pp  (when  relevant,  a  marking 
should  be  considered  a  column  vector).  Z)+,  D°,  and  g  are  “marking-dependent”, 
that  is,  they  are  specified  as  functions  of  the  marking. 

•  r  is  a  finite  set  of  transitions.  P  fl  T  =  0. 

•  Vp  £  P,Vf  €  r,V/u  £  D~t{p)  £  IN,  Dp^t{p)  £  IN,  and  Dp^^{p)  £  IN  are  the 
multiplicities  of  the  input  arc  from  p  to  t.  the  output  arc  from  t  to  p,  and  the  inhibitor 
arc  from  p  to  t,  when  the  marking  is  p,  respectively. 

•  y  C  T  X  T  IS  an  acyclic  (pre-selection)  priority  relation. 

•  ^  p  £  5't(/i)  £  {0,1}  is  the  guard  for  t  in  marking  p. 

•  £  in'^i  is  the  initial  marking. 

Places  and  transitions  are  drawm  as  circles  and  rectangles,  respectively.  The  number  of 
tokens  in  a  place  is  written  inside  the  place  itself  (default  is  zero).  Input  and  output  arcs 
have  an  arrowhead  on  their  destination,  inhibitor  arcs  have  a  small  circle.  The  multiplicity 
is  written  on  the  arc  (default  is  the  constant  1);  a  missing  arc  indicates  that  the  multiplicity 
is  the  constant  0.  The  default  value  for  guards  is  the  constant  1. 

A  transition  f  £  P  is  enabled  in  marking  p  iff  all  the  following  conditions  hold: 

1.  gt{p)  =  1. 

2.  Vp  £  P,P"t(p)  <  pp. 

3.  Vp  £  P,  D^tip)  >  pp  or  Dl  ^{p)  =  0. 

4.  ^/u  E  T.ii  t  or  n  is  not  enabled  in  p  (well  defined  because  is  acyclic). 
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Geom(a,3)  .b 

1  ■  1 ->(T>-a-^^ 

Const(2)  .  b 

1  KD-i 

2  w.p.  a 

3  w.p.  7  ■ 

=>=  w.p.  l-a-7 

1  >{3) — 1 

(l-a-7)/(l-a) 

Figure  1:  Examples  of  DDP  distributions. 

Let  S(ji)  be  the  set  of  transitions  enabled  in  marking  p.  A  transition  t  £  ^(p)  can  fire, 
causing  a  change  to  marking  M{t,  p),  obtained  from  p  by  subtracting  the  “input  bag” 
and  adding  the  “output  bag”  to  it:  =  j.L 

where  D  =  —D~  is  the  incidence  matrix.  M  can  be  extended  to  its  reflexive  and  transitive 

closure  by  considering  the  marking  reached  from  p  after  firing  a  sequence  of  transitions.  The 
reachability  set  is  given  by  1Z  =  {p  :  Bu  G  T*  A  p  =  yVf  (a,  pf®^)},  where  T*  indicates  the  set 
of  transition  sequences. 


3  Discrete  phase  distributions 

We  now  define  the  class  V  of  (possibly  defective)  discrete  phase  (DDP)  distributions,  which 
will  be  used  to  specify  the  duration  of  a  firing  time  in  a  SPN.  A  random  variable  X  is  said  to 
have  a  DDP  distribution,  X  ~  D,  iff  there  exists  an  absorbing  DTMC  ;  h  €  IN}  with 
finite  state  space  A  =  {0,1,...,??-}  and  initial  probability  distribution  given  by  [Pr{Af°^  = 
i}J  G  A],  such  that  states  A\  {0,??}  are  transient  and  states  {O.n}  are  absorbing,  and  X  is 
the  time  to  reach  state  0:  A’  =  min{A-  >  0  :  =  0}.  If  PrjA^^^  =  0}  >  0,  the  distribution 

has  a  mass  at  the  origin.  If  Pr{4f°J  =  i)  >  0  and  state  i  can  reach  state  n,  the  distribution 
is  (strictly)  defective. 

V  is  the  smallest  class  containing  the  distributions  Const(O),  Const(l),  and  Const(oo) 
and  closed  under: 

•  Finite  convolution:  if  A’l  V  and  X2  ~  P,  then  X  =  Xi  +  X2  ~  V. 

•  Finite  weighted  sum:  if  Ad  ~  P,  Ah  ~  V  and  B  G  {0,1}  is  a  Bernoulli  random 
variable,  then  X  =  BX^  +  (1  —  P)Ah  ~  P. 

•  Infinite  geometric  sum:  if  {AT  ~  P  :  A;  G  IN'*'}  is  a  family  of  iid’s  and  N  is  a  geometric 
random  variable,  then  A  =  Yli<k<N  ~  P- 
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Figure  2;  Equivalent  DTMC  representations 

The  geometric  and  modified  geometric  distributions  with  arbitrary  positive  integer  step, 
Geom(Q,co)  and  ModGeom(Q;,u;),  0  <  a  <  1,  a;  €  IN"^,  the  constant  non-negative  integer 
distribution,  Const (co),  a;  €  IN,  and  any  discrete  distribution  with  finite  non-negative  integer 
support  are  special  cases  of  DDP  distributions.  An  example  of  a  random  variable  with  non¬ 
negative  integer  support  which  does  not  have  a  DDP  distribution  is  where  N  ~  Geom(Q;). 

Fig.  1  shows  examples  of  DDP  distributions.  The  “initial  state”  5,  for  begin,  has  zero 
sojourn  time  and  is  introduced  to  represent  graphically  the  initial  probability  distribution. 
We  use  this  representation  since  it  allows  to  estimate  the  “complexity”  of  a  DTMC  by 
counting  the  number  of  nodes  and  arcs  in  its  graph.  For  simplicity,  the  last  state,  e.g., 
4  for  Geom(Q,3)  and  -3  for  Const(2),  can  be  omitted  if  it  is  not  reachable  from  b  (if  the 
distribution  is  actually  not  defective).  Unfortunately,  the  DTMC  corresponding  to  a  given 
DDP  distribution  might  not  be  unique,  even  if  the  number  of  states  if  fixed.  For  example, 
the  time  X  to  reach  state  0  for  the  DTMCs  in  Fig.  2,  both  with  five  nodes  and  seven  arcs, 
has  distribution  Unif(0,3),  that  is,  Pr{Ai  =  i]  —  1/4,  for  i  €  {0, 1,2,3}. 

4  The  DDP-SPN  formalism 

SPNs  are  obtained  when  the  time  that  must  elapse  between  the  instant  a  transition  becomes 
enabled  and  the  instant  it  can  fire,  or  firing  time,  is  a  random  variable.  By  restricting  the 
firing  times  distributions  to  P,  we  obtain  the  DDP-SPNs,  corresponding  to  a  stochastic 
process  where  the  state  has  the  form  s  =  (//,  (/>)  6  x  IN^^L  The  structural  component 
fx  is  simply  the  current  marking.  The  timing  component  (j)  describes  the  current  “phases”, 
the  state  for  the  DTMC  chosen  to  encode  the  DDP  distribution  associated  with  the  firing 
time  of  each  transition.  The  firing  time  of  a  transition  t  elapses  when  its  phase  6t  reaches  0. 
Formally,  a  DDP-SPN  is  a  tuple  T,  D~ ,  D'^ ,  D° ^  G,  F,  where: 

•  define  a  PN. 

•  Vt  G  T.'if.i  G  TZ.  C  IN  is  the  finite  set  of  possible  phases  in  which  transition  t  can 
be  when  the  marking  is  jx. 
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•  V/i  G  7?,, Vi  G  £{ii).yi,j  G  is  the  probability  that  the  phase  of  t 

changes  from  i  to  j  at  the  end  of  one  step,  when  t  is  enabled  in  marking  /i.  Hence, 

Gtil-iThj)  =  1- 

•  Vp  G  Tzyu  G  £{fi)yt  G  r,Vi  G  ^tiyyj  e  fi)),  is  the  prob¬ 

ability  that  the  phase  of  t  changes  from  i  to  j  when  u  fires  in  marking  /i.  Hence, 

1.1))  Gu, till-,  hj)  —  1- 

•  \/t  ^  T,  G  is  the  phase  of  t  at  time  0. 

•  >-  C  r  X  r  is  an  acyclic  (post-selection)  priority  relation. 

•  V//  G  TZyS  C  £{ii)yt  G  S,Wf\sili)  G  IR'*'  is  the  firing  weight  for  t  when  S  is  the  set  of 
candidates  to  lire  in  marking  /i. 

A  transition  t  G  T  is  said  to  be  a  candidate  (to  fire)  in  state  s  =  (/x,  6)  iff  all  the  following 
conditions  hold: 

1.  t  G  £{pi)- 

2.  ^t  =  0. 

3.  Vu  G  T,  n  ^  t  or  is  not  a  candidate  in  s  (remember  that  y-  is  acyclic). 

Let  C{s)  be  the  set  of  candidates  in  state  s.  Gtifi.*-*)  is  the  one-step  transition  probability 
matrix  of  the  DTMC  :  k  G  IN},  with  state  space  $t(p),  corresponding  to  the  DDP- 
distributed  firing  time  for  transition  t  in  marking  n  in  isolation,  that  is,  assuming  that  no 
other  transition  firing  affects  the  firing  time  of  t.  However,  if  another  transition  ii  fires 
before  t.  leading  to  marking  ii',  the  phase  ©;  of  t  will  change  according  to  the  distribution 
Furthermore,  after  the  firing  of  w,  the  phase  of  t  will  evolve  according  to 
Gi(^',  •,•),  which  might  differ  from  G't(//,»,«),  it  can  even  have  a  different  state  space, 
$i(/i')  instead  of  $t(^). 

We  stress  that  pre-selection  and  post-selection  have  a  different  semantic.  Only  in  the 
case  of  immediate  transitions  the  two  become  equivalent.  Assume  that  only  t  and  u  satisfy 
the  input,  inhibitor,  and  guard  conditions  in  //.  We  have  three  options,  resulting  in  three 
different  behaviors: 

•  Specify  a  pre-selection  priority  between  them,  for  example  t  y  u,  so  that  u  will  not  be 
enabled  when  t  is.  This  means  that  the  phase  <f>t  of  t  evolves  according  to  Gt(p,  •,  •), 
while  6u  does  not.  The  same  effect  would  be  achieved  using  a  guard  gu^)  =  0. 
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•  Specify  no  pre-selection  priority,  but  a  post-selection  priority  between  them,  for  exam¬ 
ple  t  'y  u.  This  means  that  the  phases  of  both  t  and  u  evolve  in  pi.  The  first  one  to 
reach  phase  0  will  fire  but,  in  case  of  a  tie,  t  will  be  chosen.  However,  if  (f)u  =  0  when 
t  fires  and  if  Ff^u(/^,0,0)  =  1,  w  might  be  a  candidate  in  the  new  marking,  and  fire 
immediately  after  t. 


•  Specify  neither  a  pre-selection  nor  a  post-selection  priority  between  them.  Then,  as  in 
the  previous  case,  t  and  u  are  in  a  race  to  reach  phase  0,  but  a  tie  is  now  resolved  by 
a  probabilistic  choice  according  to  the  the  weights:  and  «:„!{<,„}(/«),  respec¬ 

tively,  where  lu  is  a  normalization  of  w  to  ensure  that  the  weights  of  the  candidates  in 
a  marking  sum  to  one. 

Let  be  the  state  of  the  DDP-SPN  at  step  n.  Then,  the  process  :  n  G  IN} 

is  a  DTMC  with  state  space  S  C  x  Its  one-step  transition  probability  matrix 

n  is  determined  by  considering  the  possibility  of  simultaneous  firings.  Consider  a  state 
s  —  {fi.  6).  If  C{.s)  /  0,  one  of  the  candidates  will  fire  immediately,  and  the  sojourn  time  in 
s  is  zero.  Otherwise,  the  sojourn  time  in  s  is  one.  Following  GSPN  [3]  terminology,  we  call 
s  a  vanishing  or  tangible  state,  respectively.  Hence,  s  is  tangible  iff  ^  >  0. 

Let  Ss.s>  be  the  set  of  possible  event  sequences  events  leading  from  a  tangible  state 
s  =  {pi,  6)  to  a  tangible  state  s'  =  (/<',  (f>')  in  one  time  step: 


{a  =  t('), . . .  d'"))  : 

Vt  e  S{n),  Gtip,  d>u  4>t^^)  >  0, 

V2:,0  <  z  <  ^ 


(1) 

(2) 

(3) 


(1)  considers  the  one-step  evolution  of  the  phases  for  the  enabled  transitions  in  isolation, 
while  (2)  and  (3)  consider  the  sequentialized  firing  in  zero  time  of  zero  or  more  transitions 
at  the  end  of  the  one-step  period.  Hence,  (//b),  (/»W)  is  a  vanishing  state,  for  0  <  ?  <  u. 

The  value  of  the  nonzero  entries  of  H  is  obtained  by  summing  the  probability  of  all 
possible  sequences  leading  from  s  to  s': 


I]  f  n  Gt{n,<j)t,(i)T) 

■  (n  *,<oic, „(.),#.>)((■“')  fn  j 

\i=0  \teT  // 


In  a  practical  implementation,  H  is  computed  one  row  at  a  time.  The  complexity  of 
computing  row  s  of  H  can  be  substantial,  depending  on  the  length  and  number  of  sequences 
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Figure  3;  (0,  0)  can  reach  an  infinite  number  of  markings  in  one  time  step. 

in  IJs/  Ss,s’-  If  Ua'  Ss,s'  is  infinite,  special  actions  must  be  taken.  This  can  happen  for  two 
reasons: 

•  7^  is  itself  infinite,  and  state  s  can  reach  an  infinite  number  of  states  in  a  single 

step.  Consider,  for  example,  a  single  queue  with  batch  arrivals  of  size  >  0,  where 

N  ~  Geom(Q),  as  in  Fig.  3.  Following  the  firing  of  t,  a  geometrically  distributed 
number  of  tokens  will  be  placed  in  p2'-  when  the  token  is  finally  removed  from  pi  (by  the 
firing  of  v).  p2  contains  N  tokens  with  probability  a(l  — This  represents  a  batch 
arrival  of  size  N  at  the  server  modeled  by  place  p2  and  transition  y.  Unfortunately, 
finiteness  of  TZ  is  an  undecidable  question  for  the  class  of  Petri  nets  we  defined,  since 
transition  priorities  alone  make  them  Turing  equivalent  [2]. 

•  Ss.s>  can  be  infinite  for  a  particular  s'.  If  71  is  finite,  this  requires  the  presence  of 
arbitrarily  long  paths  over  a  finite  set  of  vanishing  states,  just  as  for  a  “vanishing 
loop”  in  a  GSPN  [11].  In  a  practical  implementation,  these  cycles  can  be  detected  and 
managed  appropriately. 

The  size  of  the  DTMC  underlying  a  DDP-SPN  is  affected  by  the  choice  of  the  representa¬ 
tion  for  the  DDP  distributions  involved.  Consider,  for  example,  the  DDP-SPN  in  Fig.  4(a), 
and  assume  that  transitions  U,  G,  and  have  firing  time  distributions  Const(l),  Unif(0,3), 
and  Const (2),  respectively.  The  corresponding  DTMCs  obtained  using  the  two  representa¬ 
tions  of  Fig.  2  for  Unif(0,3)  are  shown  in  Fig.  4(b)  and  4(c),  respectively.  The  number  of 
states  is  ten  in  the  first  case,  seven  in  the  second  (the  value  of  Ot  is  specified  as  whenever 
t  is  not  enabled  and  either  it  cannot  become  enabled  again  or  its  phase  is  going  to  be  reset 
upon  becoming  enabled).  The  difference  between  the  size  of  the  two  DTMCs  is  due  to  a 
lumping  [22]  of  the  states,  and  it  would  be  even  greater  if  G  had  a  more  complex  distribu¬ 
tion.  By  postponing  the  probabilistic  decision  as  much  as  possible,  the  second  DTMC  lumps 
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Figure  4:  The  effect  of  equivalent  Unif(0,3)  representations. 

states  (Oil,  •12),  (Oil,  •22),  and  (011,^32)  of  the  first  DTMC  into  a  single  one,  (011,^32), 
and  states  (011,^11)  and  (011,^21)  into  (Oil, •21). 


5  Analysis  of  DDP-SPNs 


When  using  a  SPN  to  model  a  system,  a  reward  rate  is  associated  to  each  marking  p. 
Starting  from  {(pf^^,  :  n  €.  IN},  it  is  then  possible  to  define  two  continuous-parameter 

processes:  {y{0),0  >  0},  describing  the  instantaneous  reward  rate  at  time  0:  y{0)  =  Pti(e), 
where  and  {Y{0),0  >  0},  describing  the  reward  accumulated  up  to  time 


We  consider  the  computation  of  the  expected  value  of  y{0p)  and  Y{0f)  for  finite  values 


of  6f-  Let  = 


TTi 


the  state-space  S  correspond 


|^Pr{sl”l  =  s}J  be  the  state  probability  vector  at  time  n.  Once 
ling  to  the  initial  state  (pt^l,  has  been  generated,  any  initial 
probability  vector  over  S  can  be  used  for  the  initial  probability  vector  there  is  no 
requirement  to  use  a  vector  having  a  one  in  position  and  a  zero  elsewhere.  From 

7r[°^,  we  can  obtain  tt^I  iteratively,  performing  n  matrix-vector  multiplications: 


ttI”]  =  TT^-iln 


(4) 


Since  the  DTMC  can  change  state  only  at  integer  times,  Tr{9)  —  ttLI  for  6  G  [n,  n  -f  1).  Practi¬ 
cal  implementations  assume  that  the  state  space  is  finite  and  that  the  transition  probability 
matrix  11  is  computed  before  starting  the  iterations.  The  following  shows  the  pseudo-code 
to  compute  E[y{9F)]  and  E[Y{6f)]  with  the  “power  method”: 

1.  “compute  S.  n,  and  7rl°h''; 

2.  y  4-  0;  TT  4-  ttM; 

3.  for  n  =  1  to  |_0fJ  do 
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Figure  5;  A  DDP-SPN  with  an  ergodic  underlying  DTMC. 

4.  Y  —Y  +J2{n,4,)es 

5.  It  <—  ttII; 

6.  E[Y{6f)]  <-  Y  +  {0f-  [0f\)E(h,4>)€S  Pn^{n,<t>y 

7.  E[y{0F)]  ^  J2{n,4>)es  Pfi'^{fi,<t>y 

If  the  state  space  S  is  finite,  it  is  possible  to  approximate  the  steady-state  probability  vector 
TT*  =  lim„^co  by  iterating  the  power  method  long  enough.  If  the  DTMC  is  ergodic, 
though,  other  numerical  approaches  are  preferable,  based  on  the  relation  tt*  =  tt*!!,  which 
can  be  rewritten  as  the  homogeneous  linear  system  7r*(n  —  /)  =  0,  subject  to  K  — 
Fast  iterative  methods  such  as  successive  over-relaxation  (SOR)  [12]  or  multilevel  methods 
[20]  can  then  be  employed,  although  their  convergence  is  not  guaranteed.  Fig.  5  offers  an 
example  of  an  ergodic  DTMC  obtained  from  a  DDP-SPN. 

6  Coping  with  large  state  spaces 

The  power  method  algorithm  described  requires  to  generate  the  state  space  S  and  11,  and 
then  to  iterate  using  Eq.  (4),  hence  it  assumes  a  finite  S.  However,  a  “dynamic”  state  space 
exploration  has  been  proposed  to  remove  this  restriction  [16,  18].  The  general  idea  is  to 
start  from  the  initial  state,  or  set  of  initial  states  and  iteratively  compute  both  the  set  of 
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reachable  states  and  the  probability  of  being  in  them  after  n  steps,  for  increasing  values  of  n. 
The  approach  has  been  proposed  for  the  transient  analysis  of  CTMCs  using  uniformization 
[17.  21;  27],  where,  in  practice,  the  iterations  must  be  stopped  at  a  large  but  finite  n,  thus 
resulting  in  a  truncation  error  which  can  be  bounded.  However,  the  same  approach  is  even 
more  appropriate  for  the  transient  analysis  of  DTMCs,  since,  in  this  case,  no  truncation  is 
required;  the  exact  number  of  steps, to  be  considered  is  determined  by  the  time  Op  at  which 
the  results  are  desired. 

Let  be  the  set  of  states  explored  at  step  n.  States  in  5  \  <5^”]  have  zero  probability 
at  step  n,  given  the  initial  state(s).  Then,  is  completely  determined  by  tt^,  which  is 
given,  and  <5^"^  is  obtained  from  by  considering  the  nonzero  entries  in  H*,.  for  each 

,5  ^  The  pseudo-code  for  this  modified  power  method  algorithm  is: 

1.  n  ^  0;  TT  ^  0;  0  0;  <S  ^  N  ^  5;  ^(^[o],^[o])  ^  1.0; 

2.  for  n  =  1  to  \0f\  do 

3.  Y  =  Y  +YZ(tj.,<b)£S 

4.  M'  ^  0; 

5.  while  3^  G  do 

6.  for  each  s'  such  that  Ss,s'  /  0  do 

7.  “compute  lists'”; 

8.  if  s'  ^  <S  then 

9.  A'”' ^  AT'U  {s'};  ^  5U{s'}; 

10.  A'' e- A/"\  {s}; 

11.  j\f  ^  Af'-, 

12.  TT  <—  ttH; 

13.  E[Y{0f)]  ^  Y  +  LM); 

14.  E[y{6F)]  e-  YZ(^i,d>)es  Pfi'^{M^,4>)^ 

At  the  beginning  of  the  n-th  iteration,  S  and  S\Af  contain  the  states  reachable  in  less  than 
n  and  n  -  1  steps,  respectively.  The  rows  H^^,  for  the  states  s  G  \  A/”  have  been  built  in 
previous  iterations,  while  those  corresponding  to  states  s  ^  jV  still  need  to  be  computed. 
During  the  n-th  iteration,  Af  accumulates  the  states  reachable  in  exactly  n.  but  not  fewer, 
steps.  These  states  will  be  explored  at  the  next  iteration.  This  algorithm  allows  to  study  a 
DDP-SPN  regardless  of  whether  <S  is  finite  or  not,  provided  that: 

•  6f  is  finite  (transient  analysis). 

•  A  finite  set  of  states  has  nonzero  initial  probability:  |{s  :  >  0}1  <  oo. 
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•  Each  row  of  11  contains  a  finite  number  of  nonzero  entries  or,  in  other  words,  if  the 
marking  is  fi  at  time  6.  the  set  of  possible  markings  at  time  ^  +  1  is  finite. 

The  first  two  requirements  can  be  easily  verified.  The  third  requirement  is  certainly  sat¬ 
isfied  if  Sg.s’  does  not  contain  arbitrarily  long  sequences.  This  requirement  does  not  allow 
to  analyze  exactly,  for  example,  the  DDP-SPN  in  Fig.  3.  However,  this  behavior  can  be 
approximated  arbitrarily  well  using  a  truncated  geometric  distribution  for  the  size  of  the 
batch  arrivals.  Incidentally,  we  observe  that  the  continuous  version  of  this  SPN,  where  t 
and  y  are  exponentially  distributed,  shows  that  Proposition  1  in  [9]  does  not  hold  for  un¬ 
bounded  systems:  there  is  no  SPN  with  only  exponentially  distributed  firing  times  equivalent 
to  this  GSPN  (equivalently,  there  is  no  SPN  with  only  geometrically  distributed  firing  times 
equivalent  to  the  DDP-SPN  in  Fig.  3). 

6.1  Truncating  the  state  space 

The  modified  power  method  algorithm  can,  in  principle,  perform  the  transient  analysis  of  any 
DDP-SPN  that  reaches  only  a  finite  number  of  markings  (hence  states)  in  a  finite  amount 
of  time.  In  practice,  though,  the  number  of  markings  reachable  in  a  finite  time  might  still 
be  too  large,  hence  we  need  to  find  ways  to  reduce  the  memory  requirements. 

A  first  observation  allows  us  to  reduce  the  number  of  states  that  must  be  stored  without 
introducing  any  approximation.  If  all  the  firing  times  have  geometric  distributions  with 
parameters  less  than  one,  there  is  a  nonzero  probability  of  remaining  in  a  state  a  for  an 
arbitrary  number  of  steps,  once  s  is  entered.  Indeed,  the  assumption  of  our  modified  power 
method  algorithm,  and  of  [16,  18],  is  that  the  set  of  explored  states  never  decreases:  C 

5[i]  c  c  ■  • 

However,  some  firing  times  might  have  distributions  with  finite  support,  so  it  is  possible 
that  TT^  >  0  while  =  0  and,  in  this  case,  state  s  can  be  discarded  before  computing 

Then.  we  can  redefine  <5^"^  to  be  the  set  of  time-reachable  states  at  step  n,  that  is, 
the  states  having  a  nonzero  probability  at  step  n:  =  {a:7rW  >0}. 

is  completely  determined  by  7r[°],  which  is  given,  and  hence  5^"^,  is  obtained 
from  by  computing  H^,.  for  each  s  G  and  then  restricting  the  usual  matrix- 

vector  multiplication  Tfi"]  =  to  the  entries  corresponding  to  since  the  other 

entries  are  zero  anyway.  Extreme  cases  are  illustrated  in  Fig.  6,  where,  in  (a),  = 

{(Ij.  1)  :  0  <  j  <  n},  wTile,  in  (b),  =  {(In,  1)}. 

Hence,  if  a  state  5  is  time- reachable  at  step  n,  but  time-unreachable  at  step  n  -h  1,  we 
can  destroy  it  and  its  corresponding  row  in  H  at  the  end  of  step  n  +  1.  At  worst,  the  same 
state  s  might  become  time-reachable  again  at  a  later  step,  and  the  algorithm  will  have  to 
compute  the  corresponding  row  Hj.,  for  the  transition  probability  matrix  again. 
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Figure  6:  A  case  where  C  and  another  where  ^ 

The  observation  that  the  are  not  required  to  be  a  sequence  of  nondecreasing  subsets 
is,  we  believe,  new.  Unfortunately,  geometric  distributions  with  parameter  less  than  one  are 
often  used  in  practice,  resulting  in  an  increasingly  larger  set  of  states  to  be  stored  at  each 
step  of  the  modified  power  method. 

Further  observing  that  some  markings  might  have  negligible  probability,  however,  allows 
us  to  avoid  keeping  in  its  entirety,  at  the  cost  of  an  approximate  solution,  but  with 
computable  bounds.  For  example,  in  Fig.  6(a),  the  probability  of  marking  (1/^,1)  at  step 
u,  k  <  n,  is  (^)  0.9^'0.1”“^',  which  is  extremely  small  when  the  difference  between  n  and  k 
is  large.  An  approximate  solution  approach  based  on  truncation  of  the  state-space  might 
then  be  appropriate.  At  step  n,  only  the  states  in  <5^”^  C  are  considered.  For  each  state 
s  G  its  computed  probability  is  an  approximation  of  the  exact  probability  at 
step  n: 

1.  Initially,  7r[°^  is  known,  so  set 

■(—  and  <—  {s  ;  >  0}. 

The  “total  known  probability  mass”  and  the  “total  known  sojourn  time”  at  the  begin¬ 
ning  are 

«[0]  ^  p[0]||^  =  ^  ^[0]  ^  1  and  A'M  0. 

2.  As  the  iteration  progresses,  the  size  of  might  grow  too  large  and  states  with  proba¬ 
bility  below  a  threshold  c  must  be  truncated,  destroying  them  and  their  corresponding 
row  in  11 ,  de  facto  setting  their  probability  to  zero 

for  each  s  G  (5^"^  do  if  -n-M  <  ^  tl^en  Tr^”!  <—  0. 

Compute  the  new  set  of  kept  states 

^  {s  :  ttW  >  0}. 
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Regardless  of  whether  truncation  is  performed,  the  total  known  probability  mass  at 
step  n  and  the  total  known  sojourn  time  up  to  step  n  are 


TT^ 


1  - 


E 

s€.Sl"l 


<  1 


and 


/i  W  ^  < 


n. 


Without  other  information,  we  can  only  say  that  the  probability  of  being  in  state 
5  €  at  step  n  is  at  least  ^[”1,  while  we  do  not  know  how  the  unaccounted  probability 
mass  should  be  redistributed  (we  know  that  it  should  be  redistributed  over 

the  states  in  hence  some  of  it  could  be  over  states  in  C  cSW,  but  we  have  no 
way  to  tell).  An  analogous  interpretation  holds  for 


3.  Truncation  can  be  performed  as  many  times  as  needed,  although  every  application 
reduces  the  value  of  thus  increases  our  uncertainty  about  the  state  of  the  system. 

4.  Upon  reaching  time  6p,  we  know  that,  with  probability  at  least  the  system  is 

in  one  of  the  non-truncated  states  Conversely,  a  total  of 


K{0f)  ^  Of-  -  L^fJ) 


sojourn  time  units  are  unaccounted  for.  Hence,  assuming  that  the  reward  rates  associ¬ 
ated  to  the  states  have  an  upper  and  lower  bound  pi  and  pp,  E[Y{9f)]  and  E[y[dF)] 
can  be  bounded  as  well.  If  E[Y(6f)]  and  E[y{dF)]  are  the  approximations  obtained 
using  our  truncation  approach, 

<  £[!/(«f)J  <  £|!i(P)j  +  -  fIM), 

EIY(0f)]  +  nT<i«F}  <  £|V(«f)]  <  ElY{$F)]-t-  pcKiOr)- 

Highly-reliable  systems  are  particularly  good  candidates  for  this  state-space  truncation,  since 
they  have  a  large  number  of  low-probability  states. 


6.2  Embedding  the  DTMC 

When  performing  steady-state  analysis,  it  is  possible  to  perform  an  embedding  of  the  DTMC, 
observing  it  only  when  particular  state-to-state  transitions  occur.  For  a  simple  exam¬ 
ple,  consider  the  DTMC  in  Fig.  5,  which  has  a  transition  from  state  (000111,  •  •  •!)  to 
state  (111000, 132«)  with  probability  one.  If  the  firing  time  of  transition  t4  were  changed 
to  Const(7),  instead  of  Const(l),  the  DTMC  would  have  to  contain  six  additional  states, 
(000111,  •  •  •7)  through  (000111,  •  •  *2).  This  is  obviously  undesirable,  and  it  can  be  eas¬ 
ily  avoided  by  an  embedding.  The  DTMC  of  the  embedded  process  is  exactly  that  of  Fig. 
5,  we  must  simply  set  the  expected  holding  time  of  each  state  s  to  one,  except  that 
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of  (000111,  •  •  •!),  which  is  set  to  seven.  Then,  we  can  solve  the  embedded  DTMC  for 
steady  state  and  obtain  a  steady-state  probability  vector  if  for  the  embedded  process.  The 
steady-state  probability  vector  of  the  actual  process  is  then  obtained  by  weighting  tt  accord¬ 
ing  to  the  holding  times,  a  well  known  result  applicable  to  the  steady-state  solution  of  any 
semi-Markov  process  [8]:  tt^  =  iJ2ues  ^uhu) 

For  transient  analysis,  the  same  idea  can  be  applied,  but  in  a  much  more  restricted 
way.  If,  at  step  n,  every  state  in  is  such  that  the  minimum  time  before  a  change  of 
marking  is  A-  >  1,  we  can  effectively  perform  an  embedding.  In  the  modified  power  method 
algorithm,  this  requires  advancing  n  by  k  instead  of  just  one  step  in  the  outermost  loop  and 
adjusting  the  increment  of  Y  in  statement  3  accordingly.  It  should  be  noted,  however,  that 
this  situation  is  unlikely  to  occur,  since  the  set  5^”^  may  contain  many  states  5  =  (//,  ^),  and, 
for  each  of  them,  the  DTMC  describing  the  firing  time  of  each  enabled  transition  f  in  p  must 
satisfy  min|/  G  IN  :  Pr{4^^  =  0  [  =  (ht}  >  O}  >  A:.  This  is  analogous  to  the  requirement 

for  an  efficient  application  of  adaptive  uniformization  [28]  and,  as  stated  in  the  introduction, 
it  is  unlikely  to  happen  in  general,  especially  for  large  values  of  Op. 

7  Conclusion  and  future  work 

We  defined  a  class  of  discrete-time  distributions  which,  when  used  to  specify  the  firing  time 
of  the  transitions  in  a  stochastic  Petri  net,  ensures  that  the  underlying  stochastic  process  is 
a  DTMC.  We  then  gave  conditions  under  which  the  transient  analysis  of  this  DTMC  can  be 
performed  even  if  the  state  space  is  infinite.  In  practice,  though,  the  memory  requirements 
might  still  be  excessive,  hence  we  explored  some  state-space  reduction  techniques. 

The  implementation  of  a  computer  tool  based  on  the  DDP-SPN  formalism  is  under  way. 
In  particular,  algorithms  for  the  efficient  computation  of  the  rows  of  the  transition  probability 
matrix,  II^,,,  are  being  explored. 
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